home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 22 / Cream of the Crop 22.iso / program / tn2.zip / FFT.T < prev    next >
Text File  |  1996-11-15  |  3KB  |  157 lines
  1. %
  2. % "fft.t" computes Fast Fourier Transform
  3. %
  4. %   Sample program for the T Interpreter by:
  5. %
  6. %   Stephen R. Schmitt
  7. %   962 Depot Road
  8. %   Boxborough, MA 01719
  9. %
  10.  
  11. const NU : int := 7
  12. const DIM : int := 2^NU                 % number of sample points
  13. const T : real := 0.125                 % sample interval
  14. type carray : array[DIM] of complex
  15. const PI : real := 2 * arccos( 0 )
  16. var W : carray
  17.  
  18. program
  19.  
  20.     var k, p, n : int
  21.     var f, s, t : real
  22.     var c : complex
  23.     var x : carray                      % time series
  24.  
  25.     put "time series:"
  26.     put "sec    ampl   - 0 +"
  27.     for k := 0...DIM-1 do
  28.  
  29.         s := sin( 4 * k * T )
  30.         t := k * T
  31.         put t:6:2, repeat( " ", round( 10 * ( 1 + s ) ) ) & "*"
  32.  
  33.         x[k].re := s
  34.         x[k].im := 0
  35.  
  36.     end for
  37.    
  38.     put "\ncrunching fft . . .\n"
  39.     fft( x )
  40.  
  41.     put "fft spectrum:"
  42.     put "freq   power"
  43.     p := DIM div 2
  44.     for n := p...DIM-1 do
  45.  
  46.         c.re := x[n].re
  47.         c.im := x[n].im
  48.         f := ( n - DIM ) / ( DIM * T )
  49.         put f:6:2, repeat( " ", round( cabs( c ) ) ) & "*"
  50.  
  51.     end for
  52.     for n := 0...p do
  53.  
  54.         c.re := x[n].re
  55.         c.im := x[n].im
  56.         f := n / ( DIM * T )
  57.         put f:6:2, repeat( " ", round( cabs( c ) ) ) & "*"
  58.  
  59.     end for
  60.  
  61. end program
  62.  
  63. %
  64. % fast fourier transform
  65. %
  66. procedure fft( var x : carray )
  67.  
  68.     var t : complex
  69.     var a, c, s : real
  70.     var h, i, j, k, m, n2, nu1, p : int
  71.     label another :
  72.     
  73.     n2 := DIM div 2
  74.     nu1 := NU - 1
  75.     k := 0
  76.  
  77.     for h := 1...NU do
  78.  
  79.         another:
  80.         
  81.         for i := 1...n2 do
  82.  
  83.             m := k div 2^nu1
  84.             p := bitrev( m )
  85.             a := 2 * PI * p / DIM
  86.  
  87.             c := cos( a )
  88.             s := sin( a )
  89.             j := k + n2
  90.  
  91.             t.re := x[j].re * c + x[j].im * s
  92.             t.im := x[j].im * c - x[j].re * s
  93.  
  94.             x[j].re := x[k].re - t.re
  95.             x[j].im := x[k].im - t.im
  96.  
  97.             x[k].re := x[k].re + t.re
  98.             x[k].im := x[k].im + t.im
  99.  
  100.             incr k
  101.  
  102.         end for
  103.  
  104.         k := k + n2
  105.         
  106.         if k < DIM then
  107.             goto another
  108.         end if
  109.  
  110.         k := 0
  111.         decr nu1
  112.         n2 := n2 div 2
  113.  
  114.     end for
  115.  
  116.     for k := 0...DIM-1 do
  117.  
  118.         i := bitrev( k )
  119.  
  120.         if i > k then
  121.  
  122.             t.re := x[k].re
  123.             t.im := x[k].im
  124.  
  125.             x[k].re := x[i].re
  126.             x[k].im := x[i].im
  127.  
  128.             x[i].re := t.re
  129.             x[i].im := t.im
  130.         
  131.         end if
  132.  
  133.     end for
  134.     
  135. end procedure
  136.  
  137. %
  138. % compute bit reversed index for fft
  139. %
  140. function bitrev( j : int ) : int
  141.  
  142.     var i, j1, j2, rv : int
  143.  
  144.     j1 := j
  145.     rv := 0
  146.  
  147.     for i := 1...NU do
  148.  
  149.         j2 := j1 div 2
  150.         rv := rv * 2 + ( j1 - 2 * j2 )
  151.         j1 := j2
  152.  
  153.     end for
  154.  
  155.     return rv
  156.  
  157. end function